Spin-1/2 Ji — J2 Heisenberg antiferromagnet on a square lattice: 
a plaquette renormalized tensor network study 
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We apply the plaquette renormalization scheme of tensor network states [Phys. Rev. E 83, 
056703 (2011)] to study the spin-1/2 frustrated Heisenberg J1-J2 model on an L x L square lattice 
with L=8,16 and 32. By treating tensor elements as variational parameters, we obtain the ground 
states for different J2/J1 values, and investigate staggered magnetizations, nearest-neighbor spin- 
spin correlations and plaquette order parameters. In addition to the well-known Neel order and 
coUinear order at low and high J2/ Ji , we observe a plaquette-like order at J2/ Ji ~ 0.5. A continuous 
transition between the Neel order and the plaquette-like order near ~ 0.40Ji is observed. The 
coUinear order emerges at « 0.62Ji through a first-order phase transition. 
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I. INTRODUCTION 

The search for exotic states in quantum magnets 
has been the topic of intensive research for the past 
decades. An extremely important question is when the 
conventional Neel order is destroyed, whaX kind of states 
can emerge. Frustrated antiferromagnetic spin systems, 
where the frustration from either the lattice geometry, 
or the presence of competing interactions, are candi- 
date systems to study these states. It is proposed that 
when the Neel order is destroyed by quantum fluctu- 
ations, only short-range correlations will survive, and 
the system enters a quantum paramagnetic state which 
can be described as a resonant valence bond (RVB) 
statcji The RVB state can either be a valence bond solid 
(VBS) phase, where some of the lattice symmetries are 
broken,^ or a featureless spin liquid with strong short- 
range correlations without any broken spin symmetry.'^''* 
One archetypical model to study the effect of frustra- 
tion from competing interactions is the antiferromagnetic 
(AF) J1-J2 Heisenberg model on a square lattice 
The Hamiltonian is given by. 



= Ji ^ S, 



J2 E s . 



(1) 



where Ji > and J2 > are the nearest-neighbor 
(NN) and next-nearest-neighbor (NNN) couplings, and 
the sums (ij) and {{ij)) run over NN and NNN pairs, 
respectively. Recent interests of this model have been 
revived by the discovery of Fe-based superconducting 
materialai^ where a weakened AF order can be described 
by this model with S > l/2XL-^ 

Properties of this model for S' = 1/2 in 2d have been 
studied extensively by a variety of methods, such as 
spin wave theoryj^ exact diagonalization(ED))^iii series 



expansionfi^i^"— large-iV expansion)^ functional renor- 
malization groupf^ Green's function method, pro- 
jected entangled pair states;^ etc. It is generally believed 
that in the region J2I Ji ^ 0.4, the ground state (GS) of 
the model is the Neel phase with magnetic long-range 
order (LRO). In the region J2/J1 > 0.65, spins in the 
GS are ordered at wave vector (7r,0) or (0,7r), showing 
so-called coUinear magnetic LRO. The GS in interme- 
diate region is proposed to be a quantum paramagnet 
without magnetic LRO, but the properties of this phase 
are still under intensive debate. There are several pro- 
posals for the GS, such as a columnar dimer state^^i^ 
a plaquette VBS order,§iil'^^ or a spin-liquid^-'A In the 
mean time, precise determination of the phase transition 
points is also not conclusive. Earlier series expansion 
studies^ estimate the quantum paramagnetic region is 
between 0.38 < J2/J1 < 0.62. Recent ED studyii using 
results of up to = 40 to perform finite-size extrap- 
olation estimates the transition points at J!^^ ~ 0.35Ji 
and ~ O.66J1. Meanwhile, studies by combination of 
random phase approximation and functional renormal- 
ization group find this nonmagnetic phase begins near 
J2/J1 « 0.4 - 0.45 and ends around 0.66 - 0.68^ 

Numerical studies of frustrated quantum spin systems 
present great challenges in dimensions greater than one. 
The ED method is hampered by the limitation of sys- 
tem size one can simulate. At present, the largest sys- 
tem size on the square lattice that can be simulated is 
N = 40fi^i^ Due to the minus sign problem;^ the pow- 
erful quantum Monte Carlo (QMC) method is not appli- 
cable to highly frustrated systems. In Id, the density ma- 
trix renormalization group (DMRG)^* algorithm, which 
generates matrix product states (MPS), can reach very 
high accuracy even for frustrated spin systems; however, 
direct extension of the algorithm to higher dimensions re- 
mains difficult. One promising proposal is to generalize 



the MPS to higher dimensions, the tensor network states 
(TNS))^"-2i which can serve as potential candidates for 
studying these systems. In the TNSs, the matrices are 
replaced by tensors of rank corresponding to the coordi- 
nation number of the lattice. On a 2d square lattice, the 
tensor T^-f,i{as) on site s has four indices, in addition to 
the physical index, which in the current case corresponds 
to the z-component CTs of a spin. 

Here, we should mention, according to the TNS rep- 
resentation, the rank of tensors is chosen according to 
the coordination number instead of the interaction pat- 
tern. In this way, the area law of entanglement entropy 
can be satisfied well if bond dimension D is big enough, 
especially when J2 not very large. 

Contracting over all bond indices gives the wave func- 
tion coefficient for a given spin state cti , . . . , crjv 42"^ In 
these tensor network based methods, one of the major ob- 
stacles is the computational complexity involved in the 
tensor contraction, then usually some type of approxi- 
mation is required to make the computation manageable. 
Several schemes have been proposed to facilitate the con- 
traction of the tensor networksi^"— In particular, a con- 
traction scheme based on the plaquette renormalization 
with auxiliary tensors is proposed to retain the varia- 
tional nature of the method, and it is shown that for the 
transverse Ising model, even with the smallest possible 
bond dimension (Z) = 2), non- mean- field results can be 
obtainedi^ 

In this paper, we use the TNS with the plaquette renor- 
malization scheme to study the Ji — J2 Heisenberg model 
on a square lattice. We find that even with a small bond 
dimension D = 2, it already provides a useful way to 
study the nature of the transition and estimate the value 
of the transition points. The rest of this paper is or- 
ganized as follows. In the following section, we review 
the plaquette renormalization scheme of TNS, and how 
to apply the scheme to the current model. Main results 
will be presented in Sec. IIIIl as well as some discussions. 
Sec. IIVI will give a brief summary. 

II. METHOD 

We investigate the ground state of frustrated Heisen- 
berg J1-J2 model on a square lattice, using the plaquette 
renormalized tensor network'^*. The trial wave function 
is written as 

|vi/)=^trr(rr$5Tr---)ki^2---), (2) 

where tTr indicates the tensor trace that all the tensor 
indices are summed over. Tg is rank-4 tensor on site s, 
with bond dimension D for each rank and =t or \. is 
the physical spin state. 

Explicit contraction of the tensor network is computa- 
tionally intensive. To keep the computational complex- 
ity from growing exponentially, auxiliary rank-3 tensors 
A^-j, are added to each level of the contraction process 
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FIG. 1. (a) Direct contraction of four connecting rank-4 ten- 
sors T with bond dimensions D results in a new tensor T' 
with bond dimensions ; (b) Plaquette renormalized ten- 
sor contraction via additional auxiliary rank-3 tensors A with 
bond dimensions D. The resulting tensor T' has the same 
bond dimension D as the original tensor T. 



(Fig. [T]), each transforms and truncates a pair of indices. 
A sequence of plaquette renormalizations, n = 1, 2, . . ., is 
carried out and the bond dimension of each rank is thus 
kept constant after every plaquette contraction.— In or- 
der to compute physical expectation values based on a 
TNS, one has to contract the tensors of a bra and ket 
state over their physical (e.g., spin) indices in addition 
to the bond indices of the tensors. Normally, one would 
first construct the double tensors by performing the sum 
over the physical indices, 

'^afccd ^ ^ ^i2i2fc2i2(^s)^/ijifei;i('^s)j (3) 

where the labels a, 5, c, d is a suitable combination of 
the indices of the bra (T^*) and ket (T") tensors, i.e., 
a = ii + D{i2 — 1), etc. In the calculation of the ma- 
trix element (^'|0|\I') of some operator involving one or 
several sites, similar tensors are constructed for the sites 
at which operators act weighted with a local expectation 
value {a'g jO^ |crs) . In addition, the renormalization double 
tensors can be also formed 

^abc ^ ^I2j2k2^iijiki (4) 

The bond dimension of each rank in the resulting double 
tensor becomes D = D^. This renormalization scheme 
reduces the maximum computational complexity^A to 
= D^^ for a double tensor network. 
The ground state wave function can be obtained by op- 
timizing the elements of tensors T, A for the ground state 
energy. Since the plaquette renormalization is introduced 
at the wave function level, instead of the constructed 
double tensor network, the method remains variational 
and the final energy will give a upper bound for the true 
ground state energy. We optimize the wave function us- 
ing the derivative-free Brent's method. '^^ Compared to 
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FIG. 2. (Color online) (a) The ground state energy per site as 
a function of J2/J1. The curves for L = 8 and 16 are shifted 
up by 0.05 and 0.10 for clarity; (b) The square of staggered 
magnetization as a function of J2/J1. 



FIG. 3. (Color online) (a) Extrapolated order parameters 
mo and mi as a function of J2/J1. (b) Finite-size scaling 
of M'^(tv,tv) and M^{tv,0) at J2/J1 = 0.5, where both order 
parameters mo and mi scale to zero in the thermodynamic 
limit. 



previous methods involving singular value decomposition 
(SVD)^i^ the environment of a given tensor is fully 
taken into account in the current scheme. However, the 
introduction of the renormalization A tensors at the wave 
function level effectively reduces the maximum support 
of the entanglement entropy area law in this tensor net- 
work. To reduce the number of free parameters, we im- 
pose symmetries on the trial wave function. We use a 
single plaquette, i.e. 2x2 = 4 sites as a unit cell (Fig.[T]), 
wherein tensors T on each site and auxilliary tensors Aq 
are assumed to be different. This unit is translated to 
generate a 4 x 4 unit and another set of auxilliary ten- 
sors Ai are added. This procedure is repeated until the 
full lattice is generated. Finally, the periodic boundary 
condition is appliedi^ 



III. RESULTS AND DISCUSSIONS 

We obtain the ground state wave function by varying 
the elements in the tensors T and A with D ~ 2, which 



describes a slightly entangled state beyond the product 
(mean-field) state {D — 1). Figure[5][a) shows the ground 
state energy with system sizes L — 8, 16, and 32. A 
clear cusp near J2/J1 = 0.62 is observed, signaling a 
first-order phase transition. A continuous change of the 
slope is found near J2/J1 = 0.4, probably indicating a 
continuous phase transition there. 

To study the details of the magnetic orders and the 
transition points, we compute the magnetic structure fac- 
tor, or the square of staggered magnetization at wave 
vector q, defined as 

^^'(^) = ]7^E«^'''^""'"'"^'^(S,.S,), (5) 

where Vi = {xi,yi), and q = (tt, tt) for the Neel order, and 
(0, 7r) or (77,0) for the collinear order. M^(q) tends to 
the square of the order parameter in the thermodynamic 
limit if there is magnetic ordering at wave vector q, and 
scales like 1/A^ in a magnetically disordered phase. 

Figure ^h) shows the results of the square of stag- 
gered magnetizations M^(7r,7r) and M^(7r,0). From the 
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FIG. 4. (Color online) The ^(black) and xyijeA) components 
of the square of staggered magnetization and the sum of the 
two (green) as a function of J2/J1. (Inset) Same quantities 
in the regime of J2/J1 = 0.45 ~ 0.65. 




FIG. 5. (Color online) The NN spin-spin correlations (Si ■ Sj) 
(black numbers near bond) and the plaquette order parameter 
(red numbers in italic) for J2/J1 = (a) 0.10 and (b) 0.50, with 
system size L = 32. We show only one corner (4 x 4) of the 
entire lattice as the pattern is repeated periodically. 



small J2/J1 side, the Neel order is smoothly suppressed 
as J2 increases, until J2/J1 — 0.40, where a discontin- 
uous jump of the Neel order is observed for L = 8, 
and the jumps become less pronounced as the system 
size increases. This strong size dependence of the jump 
is another example that in a finite-size tensor network 
state with finite bond dimensions, there exists two en- 
ergy minima near the transition, rendering the transition 
first-order at small N . For a putative continuous tran- 
sition, these two minima move closer to each other with 
increasing N and the transition becomes continuous at 
N ^00^ 

From the large J2 / Ji side, the collinear order also de- 
creases smoothly, until J2/J1 — 0.6 where a clear first- 
order transition occurs. Unlike the previous case, the 
jumps in M^(7r, 0) remain robust upon increasing TV, 
strongly suggesting against a continuous transition here. 
This transition to the collinear order is consistent with 
previous numerical calculations! "'^^i"'^*^'^"" — 



We now use our data from different sizes to ex- 
tract the order parameters in the thermodynamic limit. 
This allows us to estimate the transition points between 
the Neel/coUinear state and the non- magnetic (disor- 
dered) phase. The finite-size extrapolation rules for the 
two-dimensional antiferromagnetic Heisenberg model are 
well-known^^Ti^. Following Refs. |l3 and |4l| . we define 
the Neel order parameter as toq = 2 liniAr^oo M(7r, tt). 
This normalization is chosen so that mo = 1 in a perfect 
Neel state. The finite-size behavior of M^(7r,7r) is given 



0.62075 c 
pL 



(6) 



where c is the spin-wave velocity and p is the spin stiff- 
ness. The order parameter for the collinear order is de- 
fined as nil — \/8 hmAT^oo Af(7r, 0). The finite-size be- 
havior of M (tt, 0) is given by^^i^ 

„, , 1 o const. , , 

M2(^,0) = -m? + — ^ + ... (7) 

The extra 1/2 factor comes from the fact that the ground 
state has an extra two-fold degeneracy q — (tt, 0), (0, tt), 
and this symmetry is broken in the thermodynamic limit. 
Figure [Ha) shows the extrapolated results for mo and 
mi as a function of Jij J\. We find that the GS near 
J2/J1 = 0.5 is magnetically disordered, i.e., both mo 
and mi vanish. Figure |31[b) shows the finite-size scal- 
ing of M2(7r,7r) and M'^{t:,Q) at J2/J1 = 0.5, which 
both shows a \/N scaling with the zero intercept as 
— ?► 00. The transition points are estimated to be 
= O.4OJ1 and J^^ = 0.62Ji, consistent with esti- 
mates from series expansioni^ii^"— 1^ where J!^^ « 0.38 Ji 
and J 2^ « 0.62Ji, and slightly different from ED re- 



sults 



O.35J1 and 



0.66 Ji. ^4 Near 



we 



fit the Neel order parameter mo to a power law mo ~ 
(J2 — J^^)^, and an asymptotic mean-field behavior con- 
sistent with (3 = 1/2 is also observed.'^® For J2 = 0, we ob- 
tain mo = 0.592 which is slightly lower than the best es- 
timate from the quantum Monte Carlo (mo — 0.6140)4^ 
Although it is also possible to extract c and p from our 
data based on Eq. ([6]), it is argued that determination of 
these quantities by fitting the prefactors of the leading 
finite-size corrections (0(1/L)) can not reach the same 
accuracy as the magnetic order parameters-^"* 

Analogous to how mean-field theory produces 
symmetry-broken states, this method can produce so- 
lutions which break spin-rotation symmetry on a finite 
lattice."^*'^"^ We examine the spin-rotation symmetry of 
the ground state, with the focus in the nonmagnetic 
phase. Figure [5] shows z and xy components of the square 
of staggered magnetization at g = (tt, tt) for L = 32, de- 
fined as 



1 



T[{xi-Xj) + (yi-yj)] ^5^2 ^ 



M. 



xy 



iTT[{xi-Xj) + (yi-yj)] 
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For reference, the sum of the two is also included. In 
the Neel phase, the spin-rotational symmetry is clearly 
broken.— Increasing J2 through a phase transition to 
the strongly frustrated regime (i.e., 0.45 < J2/J1 
0.60), the spin-rotation symmetry is restored with = 
\Mly = iM^, as expected. 

In order to clarify the possible new phase in the highly 
frustrated region around J2/J1 = 0.5, we calculate the 
nearest- neighbor spin-spin correlations for L = 32. Fig- 
ure [S] shows the results for J2/J1 = 0.10, which is deep 
inside the Neel phase, and J2/J1 = 0.50, which is in the 
magnetically disordered phase. The numbers in black 
near the bond are the NN spin-spin correlation, and the 
thickness of the bond is proportional to its magnitude. 
For J2/J1 = 0.50 [Fig. [5jb)], the NN spin-spin correla- 
tions within a single plaquette are much stronger than 
those between plaquettes. On the other hand, deep in- 
side the Neel phase J2/J1 ^ 0.10 [Fig. Elja)], the NN 
spin-spin correlations shows a more uniform pattern, al- 
though weaker correlations are present in some bonds be- 
tween plaquettes. Overall, it is clear that the correlations 
inside a 2 x 2 plaquette become stronger upon increasing 
J2 / Ji , which indicates a possible plaquette order in the 
magnetically disordered phase. 

We also investigate the plaquette order parameter, 
which distinguishes clearly a Neel ordered phase from 
a plaquette order, defined as-- 

Qal3-iS = \{Pal3'yS + -Pq^V^ ^ 2[(Sq • S^)(S-y • S^) 
+ (Sa • S5)(S^ • S^) — (Sq • S^)(S^ • 85)] 

+ i(S„-S^ + S^.S5 + i). (8) 

The results of the plaquette order parameter are shown 
also in Figs. [S] (numbers in red italic) for J2/J1 — 0.10 
and 0.50. In the most frustrated region, we observe sig- 
nature of the plaquette order. For J2/J1 = 0.50, the 
plaquette order parameter is much stronger within a pla- 
quette, consistent with observation from the spin-spin 
correlations. This order parameter is small in Neel phase 
(J2/J1 — 0.10), although some traces of the plaquette 
order is still present. This might be due to the inher- 
ent structure of the renormalization scheme, which ex- 
plicitly breaks the translational invariance, or possibly 
the plaquette correlations already start to build up in 
this regime. It remains to further explore whether this 
plaquette order is favored due to our renormalization 
scheme. The plaquette renormalization scheme reduces 



the amount of entanglement support between plaquettes 
by a factor of D compared with the exact contraction. 
This may bias toward those correlations compatible with 
the plaquette structure. 



IV. CONCLUSION 

We use the plaquette renormalization scheme to study 
spin-1/2 frustrated Heisenberg J1-J2 model on a square 
lattice with different sizes of L = 8, 16, and 32. Using the 
smallest possible bond dimension D = 2 for the underly- 
ing tensors, we are already able to obtain results beyond 
the mean-field theory. Since our method is variational, 
and the calculations are done on finite lattices, we are 
able to perform finite-size scaling to extrapolate the or- 
der parameters in the thermodynamic limit. We observe 
signatures of a continuous transition at ~ 0.40 Ji, and 
a first-order phase transition at J^^ ~ 0.62Ji, consistent 
with previous numerical calculations . ^^1^"'^ Our calcula- 
tions on the NN spin-spin correlation and the plaque- 
tte order parameter indicates a possible plaquette VBS 
order for < J2 < ■ The effects of the plaque- 
tte renormalization scheme and the bond dimension D 
dependence of the physical observables require further 
studies and will be presented in a future work4^i^ 
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note added.- After submitting this manuscript, we re- 
cently learned of the DMRG work by Jiang et al^ and 
the tensor product state approach by Wang et al^ on 
the same model, which argue that the ground state in 
the nonmagnetic regime near J2/J1 ~ 0.5 could be a Z2 
spin liquid. 
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